#include "math.h"
#include "mylib.h"
#include <iostream>
using namespace mylib;
double lin_interp(double x,double x0,double x1,double y0,double y1){
    return  y0 + (y1-y0)/(x1-x0)*(x-x0);
}

double aitken(double x, VecDoub_I xi, VecDoub_I yi){
    int m;
    VecDoub tmp=yi;
    m = xi.size();


    for (int i = 1; i < m; i++)
    {
        for (int j = 0; j < m-i; j++)
        {
           tmp[j]= ((x-xi[j])*tmp[j+1] - (x-xi[j+i]) * tmp[j])/(xi[j] - xi[j+i]); 
        }
        
    }
    return tmp[0];
}

int main(int argc, char const *argv[])
{
    double x[10];
    double y[10];
    double xx;
    double yy;

    //生成插值的数据
    for (int i = 0; i < 10; i++)
    {
        x[i] = 2 * 3.1415926 * i / 10;
        y[i] = sin(x[i]);
    }

    std::cout<<"线性插值"<<std::endl;
    for (int i = 0; i < 9; i++)
    {
        xx = (x[i] + x[i + 1]) / 2;
        yy = lin_interp(xx, x[i], x[i + 1],y[i],y[i+1]);
        std::cout << "x = " << xx << ",插值y(x)=" << yy << ",Sin(x)="<<sin(xx)<<std::endl;
    }

    std::cout<<"三点二次插值"<<std::endl;
    VecDoub tmpx(3),tmpy(3);
    for (int i = 0; i < 9; i++)
    {
        if (i<=6)
        {
            for (int k = 0; k < 3; k++)
            {
                tmpx[k] = x[i+k];
                tmpy[k] = y[i+k];
            }
            
        }
        else{
            for (int k = 0; k < 3; k++)
            {
                tmpx[k] = x[6+k];
                tmpy[k] = y[6+k];
            }
        }
        
        xx = (x[i] + x[i + 1]) / 2;
        yy = aitken(xx, tmpx,tmpy );
        std::cout << "x = " << xx << ",插值y(x)=" << yy << ",Sin(x)="<<sin(xx)<<std::endl;
    }


    return 0;
}
